Snowflakeで特定エリア内の地理空間オブジェクトをカウントしてみた

モダンデータスタック(MDS)

#Snowflake

大高大輔

2021.04.15

この記事は公開されてから1年以上経過しています。情報が古い可能性がありますので、ご注意ください。

こんにちは！DA(データアナリティクス)事業本部インテグレーション部の大高です。

Snowflakeでは、2021年4月現在プレビュー機能ではありますが「地理空間データ型」および「地理空間関数」が利用できます。

「地理空間関数」を利用すると「GEOGRAPHY」型の地理空間オブジェクトに対して、重なりの検出や距離の測定などができます。今回はこの「地理空間関数」の「ST_CONTAINS」関数を利用して「特定エリア内において、地理空間オブジェクトがいくつあるのか」をカウントしてみたいと思います。

想定シナリオ

今回、データとして利用するものは以下になります。

日本の「各都道府県」のポリゴンデータ
日本全域にわたる地理空間オブジェクト

これを利用して「各都道府県内に、カウント対象の地理空間オブジェクトがいくつ存在するか」をカウントします。

利用するデータ

今回利用するデータはそれぞれ以下とします。

日本の「各都道府県」のポリゴンデータ

以下で公開されているファイルの「japan.geojson」を利用します。このファイルを事前にダウンロードしておきます。

なお、こちらのデータの出典元は「地球地図日本」のデータとなり、こちらのデータをSnowflakeにロードして利用します。

地球地図日本 | 国土地理院

日本全域にわたる地理空間オブジェクト

以下で公開されているOokla社によるネットワークスピードテストのデータを利用します。元データは日本を含めて世界全域をカバーしているので、このデータで要件を満たしているものとします。

このデータは、以下のエントリで利用したデータと同じものとなります。

前提条件

Snowflakeのウェアハウスは「Large」サイズのウェアハウスを利用します。

テストデータの準備

では、まずはテストデータを準備していきます。

都道府県ポリゴンのデータ

ダウンロードした「japan.geojson」ファイルをSnowSQLを利用して、ユーザーステージにアップロードします。

$ snowsql
❄ SnowSQL ❄ v1.2.14
Type SQL statements or !help
foo_bar#(no warehouse)@(no database).(no schema)>PUT file://./japan.geojson @~/spatial/;
japan.geojson_c.gz(1.18MB): [##########] 100.00% Done (0.175s, 6.78MB/s).       
╒═══════════════╤══════════════════╤═════════════╤═════════════╤════════════════════╤════════════════════╤══════════╤═════════╕
│ source        │ target           │ source_size │ target_size │ source_compression │ target_compression │ status   │ message │
╞═══════════════╪══════════════════╪═════════════╪═════════════╪════════════════════╪════════════════════╪══════════╪═════════╡
│ japan.geojson │ japan.geojson.gz │    13050092 │     1240891 │ NONE               │ GZIP               │ UPLOADED │         │
╘═══════════════╧══════════════════╧═════════════╧═════════════╧════════════════════╧════════════════════╧══════════╧═════════╛
1 Row(s) produced. Time Elapsed: 1.215s

その後、テーブルを作成してロードします。ロードする際ですが、元データがgeojson形式なので一旦VARIANT型のテーブルにCOPYして、その後にLATERAL FLATTENを利用してネストの解除をしつつ、実データ用テーブルにINSERTしています。

-- 都道府県のエリアテーブル（COPY用）を作成
CREATE OR REPLACE TABLE japan_geojson_raw (
  src VARIANT
);

-- 都道府県のエリアテーブル（実データ用）を作成
CREATE OR REPLACE TABLE japan_geojson (
  code NUMBER,
  name VARCHAR,
  area GEOGRAPHY
);

-- データをユーザーステージからCOPY
COPY INTO
  japan_geojson_raw
FROM
  @~/spatial/
FILE_FORMAT = (TYPE = 'JSON')
;

-- ネストを解除しつつ、実データテーブルにINSERT
INSERT INTO japan_geojson
SELECT
  value:properties.id::number AS code,
  value:properties.nam_ja::varchar AS name,
  TO_GEOGRAPHY(value:geometry) AS area
FROM
    japan_geojson_raw
  , LATERAL FLATTEN(INPUT => src:features)
;

カウント対象の地理空間オブジェクト（タイルデータ）

OOLKAのデータはS3でパブリックに公開されているので、外部ステージを作成してデータをCOPYします。

-- ステージを作成
CREATE OR REPLACE STAGE ookla url = 's3://ookla-open-data/';

-- 登録先テーブルを作成
CREATE OR REPLACE TABLE oolka_geo (
  avg_d_kbps INTEGER,
  avg_u_kpbs INTEGER,
  avg_lat_ms INTEGER,
  tests      INTEGER,
  devices    INTEGER,
  quadkey    VARCHAR,
  tile       GEOGRAPHY
);

-- parquetデータをCOPY
COPY INTO oolka_geo FROM (
  SELECT $1:avg_d_kbps::INTEGER,
         $1:avg_u_kpbs::INTEGER,
         $1:avg_lat_ms::INTEGER,
         $1:tests::INTEGER,
         $1:devices::INTEGER,
         $1:quadkey::VARCHAR,
         ST_GEOGRAPHYFROMWKT($1:tile::VARCHAR)
    FROM @ookla/parquet/performance/type=mobile/year=2020/quarter=4/2020-10-01_performance_mobile_tiles.parquet
)
FILE_FORMAT = (TYPE=PARQUET COMPRESSION=SNAPPY);

これで、データの準備は完了です！

特定エリア内の地理空間オブジェクトをカウントしてみる

では、ST_CONTAINS関数を利用してカウントしてみましょう。各都道府県上にOoklaの地理空間オブジェクトがいくつ存在するかをカウントしてみます。

なお、元データはそれぞれ以下のような件数となります。

-- 都道府県は47件
SELECT COUNT(area) AS rec_count FROM japan_geojson;
-- 47

-- Ooklaの地理空間オブジェクト数は約420万件
SELECT COUNT(tile) AS rec_count FROM ookla_geo;
-- 4,197,854

単純にSQLでカウントする

まずは愚直にSQLを書いてカウントしてみます。（もっと良い書き方があるかもしれません）

また、県コードでORDER BYを入れたほうが結果が分かりやすいのですが、実行時間測定の影響を考慮して省いています。

-- エリア内のタイルをカウント
SELECT
  japan_geojson.code,
  japan_geojson.name,
  COUNT(ookla_geo.tile) AS contains_counts
FROM
  japan_geojson, ookla_geo
WHERE
  ST_CONTAINS(japan_geojson.area, ookla_geo.tile)
GROUP BY
  japan_geojson.code,
  japan_geojson.name
;

カウント結果は以下のようになり、クエリ実行時間は「5分21秒」でした。321秒なので、1都道府県あたり6.8秒ぐらいですね。

また、カウント結果自体は都道府県の面積から考えても、合っていそうです。

CODE	NAME	CONTAINS_COUNT
3	岩手県	644
45	宮崎県	358
18	福井県	354
14	神奈川県	2930
11	埼玉県	3085
4	宮城県	1070
21	岐阜県	902
46	鹿児島県	473
43	熊本県	581
23	愛知県	3489
17	石川県	530
26	京都府	862
27	大阪府	2588
30	和歌山県	273
32	島根県	244
42	長崎県	328
44	大分県	446
5	秋田県	425
33	岡山県	798
8	茨城県	1570
39	高知県	204
20	長野県	1025
40	福岡県	1516
28	兵庫県	1828
6	山形県	535
22	静岡県	1891
24	三重県	794
47	沖縄県	530
10	群馬県	1110
13	東京都	3537
36	徳島県	264
25	滋賀県	587
41	佐賀県	254
2	青森県	627
16	富山県	499
34	広島県	988
19	山梨県	460
29	奈良県	606
9	栃木県	1070
1	北海道	2135
12	千葉県	2578
38	愛媛県	531
35	山口県	546
7	福島県	811
37	香川県	361
15	新潟県	975
31	鳥取県	267

クエリプロファイルを確認してみます。想定通りですが、JOINの方法がCartesianJoin(CROSS JOIN, デカルト積)となっておりコスト掛かっていることがわかります。

今回は対象エリアが都道府県だったので47件で済みましたが、これがN件になると更に時間がかかりそうです。

ストアドプロシージャを利用してカウントする

次に、ストアドプロシージャを利用してカウントしてみます。ストアドプロシージャでは、都道府県レコードでループを回し、各ループ内で都道府県毎にカウントして、結果テーブルにINSERTするという手法をとってみます。

1レコード毎にINSERTするのは悪手ではありますが、今回はループ件数も少ないのでこれで試してみます。ループ件数が多い場合は一括でINSERTできるようにしたほうが良いですね。

-- 結果登録テーブルを作成
CREATE OR REPLACE TABLE tile_count_result (
  code NUMBER,
  name VARCHAR,
  count NUMBER
);

-- カウント用プロシージャを作成
CREATE OR REPLACE PROCEDURE update_tile_count()
  RETURNS FLOAT NOT NULL
  LANGUAGE javascript
  AS
  $$
    // クエリ
    const select_pref_query = 'SELECT code, name, ST_ASWKT(area) AS area FROM japan_geojson;';
    const insert_result_query = 'INSERT INTO tile_count_result (code, name, count) VALUES (?, ?, ?);'
    const select_count_query = 'SELECT COUNT(ookla_geo.tile) AS count FROM ookla_geo WHERE ST_CONTAINS(TO_GEOGRAPHY(?), ookla_geo.tile);'

    // 県でループ
    const select_pref_statement = snowflake.createStatement( {sqlText: select_pref_query} );
    const select_pref_result_set = select_pref_statement.execute();
    while (select_pref_result_set.next()) {
      const pref_code = select_pref_result_set.getColumnValue('CODE');
      const pref_name = select_pref_result_set.getColumnValue('NAME');
      const pref_area = select_pref_result_set.getColumnValue('AREA');

      // 県域内のタイルをカウント
      const select_count_statement = snowflake.createStatement(
                                        {
                                          sqlText: select_count_query,
                                          binds:[pref_area]
                                        }
                                      );
      const select_count_result_set = select_count_statement.execute();
      while (select_count_result_set.next()) {
        const contains_count = select_count_result_set.getColumnValue('COUNT');

        // 結果テーブルにインサート
        const insert_result_statement = snowflake.createStatement(
                                          {
                                            sqlText: insert_result_query,
                                            binds:[pref_code, pref_name, contains_count]
                                          }
                                        );
        insert_result_statement.execute();
      }
    }

    return 0;
  $$
;

-- プロシージャを実行
CALL update_tile_count();

プロシージャ実行時間は「4分8秒」でした。少し改善しましたね！248秒なので、1都道府県あたり5.3秒ぐらいですね。INSERT処理の部分を一括でINSERTできるように修正すれば、もう少し早くなるかもしれません。

クエリプロファイルについては、ストアドプロシージャなので詳細は確認できませんが、おそらくCROSS JOINではなくなった分、パフォーマンス改善したのではないかと考えられます。

なお、カウント結果自体は以下のとおりで単純にクエリ実行したものと同じです。

-- カウント結果を確認
SELECT
  code,
  name,
  count
FROM
  tile_count_result
ORDER BY
  code
;

CODE	NAME	COUNT
1	北海道	2135
2	青森県	627
3	岩手県	644
4	宮城県	1070
5	秋田県	425
6	山形県	535
7	福島県	811
8	茨城県	1570
9	栃木県	1070
10	群馬県	1110
11	埼玉県	3085
12	千葉県	2578
13	東京都	3537
14	神奈川県	2930
15	新潟県	975
16	富山県	499
17	石川県	530
18	福井県	354
19	山梨県	460
20	長野県	1025
21	岐阜県	902
22	静岡県	1891
23	愛知県	3489
24	三重県	794
25	滋賀県	587
26	京都府	862
27	大阪府	2588
28	兵庫県	1828
29	奈良県	606
30	和歌山県	273
31	鳥取県	267
32	島根県	244
33	岡山県	798
34	広島県	988
35	山口県	546
36	徳島県	264
37	香川県	361
38	愛媛県	531
39	高知県	204
40	福岡県	1516
41	佐賀県	254
42	長崎県	328
43	熊本県	581
44	大分県	446
45	宮崎県	358
46	鹿児島県	473
47	沖縄県	530